Successes and Failures of Kadanoff-Baym Dynamics in Hubbard Nanoclusters 



M. Puig von Friesen, C. Verdozzi, and C.-O. Almbladh 

Mathematical Physics and European Theoretical Spectroscopy Facility (ETSF), Lund University, 22100 Lund, Sweden 

(Dated: May 13, 2009) 

We study the non-equilibrium dynamics of small, strongly correlated clusters, described by a 
Hubbard Hamiltonian, by propagating in time the Kadanoff-Baym equations within the Hartree- 
Fock, 2^^ Born, GW and T-matrix approximations. We compare the results to exact numerical 
solutions. We find that the T-matrix is overall superior to the other approximations, and is in good 
agreement with the exact results in the low-density regime. In the long time limit, the many-body 
approximations attain an unphysical steady state which we attribute to the implicit inclusion of 
infinite order diagrams in a few-body system. 
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Introduced almost fifty years ago, the Kadanoff-Baym 
equations (KBE) [1, 2 are a cornerstone of a microscopic 
theory of non-equilibrium quantum processes [3l [4] . In 
recent years, the KBE have received growing attention 
by the strongly correlated systems (SCS) and quantum 
transport (QT) communities, due to increased computa- 
tional capabilities [SlISllTlElElllOlIIIlIIlIIS]. Non- 
perturbative treatments of SCS with the KBE are be- 
ginning to appear, albeit specialized to spatially uniform 
fields [13 . At present, for general time-dependent (TD) 
and inhomogeneous fields, the only available treatments 
of the KBE for SCS rely on many-body perturbation the- 
ory (MBPT), where one systematically builds approxi- 
mations for the one-particle Greens function, G. The 
latter is the key quantity in the KBE. Most work with 
the KBE has been so far for model systems but ah initio 
studies based on MBPT+KBE are beginning to appear 

laiioiin]. 

Since MBPT+KBE currently appears to be the only 
route to SCS for general perturbations, it is rather sur- 
prising that a rigorous assessment of the applicability 
of the approach has not yet been made. In this work, 
we perform such an evaluation, by comparing different 
MBPT schemes for the KBE against exact solutions for 
small, isolated clusters exposed to TD inhomogeneous 
perturbations. Clearly, in finite closed systems, one can- 
not address important issues such as the onset of a steady 
state regime, or the role of dissipation: This requires cou- 
pling the clusters to an environment as done, for example, 
in [To] for the stationary limit, or in [12] for the TD case. 
However, as shown here, studying small isolated clusters 
can reveal crucial and unexpected finite size effects in the 
nature of the KBE solutions. 

Our clusters are small ID chains with Hubbard inter- 
actions. For the MBPT+KBE, we considered four con- 
serving [1] many-body approximations (MB As): Hartree- 
Fock, 2^^ Born, GW [14] and T-matrix [I5l [16] (HFA, 
BA, GWA and TMA, respectively). We have also used a 
spin dependent GWA (SGWA) [TDl [H] to mitigate spu- 
rious self-interaction effects [HI [19] . 

Another common approach to microscopic TD phe- 



nomena is Time Dependent Density Functional Theory 
(TDDFT) [20]. A well established connection exists be- 
tween TDDFT and MBPT on the Keldysh contour [21]. 
Here, to highlight this link, we will obtain MBPT-based 
exchange-correlation (xc) potentials via TD reverse engi- 
neering [22], using the TD densities from the KBE. 

The main findings of this paper are (i) for Hubbard- 
like interactions, the TMA performs very well at low 
densities, both for ground-state and time dynamics, and 
is in general superior to the other MBAs at all elec- 
tron densities. The difference in performance among 
the MBAs increases at larger interaction strength U (ii) 
the SGWA performs better than the GWA; however, the 
SGWA breaks down when U exceeds a critical value iii) 
a comparison of exact vs MBPT-based xc potentials for 
TDDFT shows trends similar to those for the densities 
iv) during the time evolution, the KBE develop an un- 
physical steady state solution. This is a central result of 
the present work, and we argue that this problem occurs 
in general, whenever MBPT is applied to finite systems, 
and self energies based upon infinite partial summations 
are used. 

Model system. The Hamiltonian of our open-ended Hub- 
bard chains is (we set the onsite energies equal to 0) 

H= -V^a^^^aR'cT^U ^fiR^fiRi +^ wr (t) fiRg (1) 

{RR')a R Ro- 

Here uro- = a\^^aRcj^ (J =T, j, and {RR') denotes nearest 
neighbor sites. The hopping parameter V = 1 and w (t) 
is a local external field which can be of any shape in time 
t and space. U and w (t) are given in units of V. We 
consider chains of lengths L = 2, 6 and = 2,6 number 
of electrons; we take spin-up and -down electrons equal 
in number, = N^; this holds at all times, since H has 
no spin-flip terms. Henceforth, n = /L. 
Kadanoff-Baym equations. We find the non-equilibrium 
one-particle Greens function 6(^1,^2) by solving 
its equation of motion {idt^ — h {ti)) G {ti^t2) = 
S (^1, ^) G (t, t2) dt [U I2T], and the one corresponding 
to ^2- Here h is the single particle Hamiltonian, S is the 
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FIG. 1: Contours and approximate self-energies 



self energy and 7 represents the Keldysh contour . In- 
stead of using the more conventional contour of Fig. [l} 
we use an equivalent one (Fig. [iji) [4 , numerically more 
stable and with an analytical limit when the temperature 
T = /^-^ ^ 0. The KBE for ti then become 

{idt,-h{h))G^ {hM) = 
Jo ^ 



1 

^ Jo 

- / dr^> {ti,ir)G< {ir,t2) 
^ Jo 



(2) 



Here ^ and R/A have the usual meaning [l|i2j. We work 
at T = 0. The initial state is the correlated [3] ground 
state, obtained by solving the Dyson equation G = Go + 
GoS[G]G self consistently, with (e - h)Go = 1. Eq. ^ 
is solved by time propagation using a predictor-corrector 
method similar to that presented in [23 . Conservation 
laws and time reversal symmetry are obeyed up to arbi- 
trary accuracy. In a finite system, exact and approximate 
spectral functions are meromorphic. For example, for the 
Greens function, Grr^ (e) = ~ where 

aj are pole positions and ^j^^/ are residue matrices in 
the single particle orbital representation. This represen- 
tation is very convenient as convolutions become simple 
matrix products [25 . 

Many-body approximations (MB As). In the Hub- 
bard model, the interaction can be treated either as 
spin-dependent, U^^nR^riRi or as spin-independent, 

\^HRcTa' ^^Rct^^Rct'^R^'^R^' Thcsc two ways are evi- 
dently equivalent in any order by order expansion (Fig. 
[1]) such as the HFA or BA. In approximations, how- 
ever, this equivalence may be lost. For example, in 
spin-independent GWA, S (12) = G(12)W(12), where 
W = U ^UPW and P(12) = G(12)G(21). In spin- 
dependent GW (SGWA), instead, W = UPU^{UPf W. 
The TMA, where S (12) = /J^G (21) T (12), is treated 
only as spin-dependent, with T = (j) — U (j)T and (12) = 
G(12)G(12). 

Exchange- correlation potential. From the TD densities 
we obtained via reverse engineering the corresponding 
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FIG. 2: (Color online) Ground state spectral functions at 
L—6 for different fillings and interactions. The curves corre- 
spond to exact (black), TMA (red), BA (green), GWA (blue) 
and HFA (brown). The curves are shifted for clearer compar- 
ison and we have used a Lorentzian broadening F = 0.2V. 



vks = vh^v^^c, vh and 

Vxc being the Hartree and the xc 
potential. In practice we minimized J dt \n {t) — uks (^)| 
[22] , where uks is the Kohn-Sham density, found by solv- 
ing i'^KS = {t ^ w ^ Vks) i^KSi and i is the kinetic term. 

Ground state. We start by solving the Dyson equation 
self consistently to obtain the initial G. Fig. [2] shows se- 
lected spectral densities for different particle concentra- 
tions and interactions. The case of half filling is shown 
in panels 2a) and 2b), where all the curves display the 
electron-hole symmetry. The exact solution exhibits the 
opening of a correlation gap on increasing U ; this is also 
reproduced, at different extent, by the approximate treat- 
ments. The gap increase is most closely reproduced (but 
largely underestimated) by the BA, with the GWA and 
the TMA giving even smaller gaps. The exact lower and 
upper Hubbard bands are incorrectly reproduced by all 
the MBAs. In particular, for [7=4, the MBAs intro- 
duce spurious satellite structures (less pronounced in the 
TMA) away from the band region. Results for the low 
density regime are in panels 2c) and 2d). For t/=l, the 
main effect of the interaction is an increased asymmetry 
in the band region; the overall agreement between exact 
and MBAs is rather good, especially for the BA and the 
TMA. For [7=4, the most notable feature in the exact 
solution is a satellite at about 6.5 (in an extended sys- 
tem, this would be a two-electron anti-bound state out- 
side the continuum). Such satellite is well reproduced by 
the TMA, smeared out in the BA and the GWA, and ob- 
viously absent in the HFA. In the band region, for U =4, 
the agreement is moderate. We find that the spectral 
functions are better at first iteration. This corresponds 
to the known fact that self-consistency often deteriorates 
the one-electron spectral properties. Self-consistency is 
required for dynamically conserving approximations and 
for total energies, but other summation criteria should 
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FIG. 3: (Color online) TD densities (panels a to e) and vks 
(panel f) on site 1. Exact (thick solid black), TMA (dashed 
red), BA (dotted green), GWA (thin solid blue) and HFA 
(brown dashed dot) results. 



be adopted for spectral densities [2T, and vertex correc- 
tions are often required to remove self-consistency arti- 
facts [25]. We found that the SGWA is slightly better 
than the GWA, since it includes fewer faulty diagrams. 
Yet, it is still worse than the BA or the TMA. It is worth 
noting that the SGWA has a magnetic instability as a 
function of U. In the dimer, where the poles of VK[Go] are 
e = ±V4F2 ± 2VU, this occurs for U > 2V. As a conclu- 
sive remark, self-consistent partial sums imply infinitely 
many diagrams, thereby violating the particle number 
constraint for finite systems. This results in an infinite, 
but discrete, number of poles in the spectral function. 
Even if this has no major effect in the ground state, it 
can have rather startling consequences for the TD behav- 
ior of the clusters. 

Time dependence T. Densities. The system is in the 
ground state for t < 0. At t > an external field w 
is applied. All results in the paper were obtained with 
a step-like w (t) = woQ{t) (time unit = |V^|~^) applied 
only to the first {R = 1) site, but we also examined other 
space and time dependencies. Fig. [3] shows the TD den- 
sities and Vks at = 1 for some representative cases. 
In panel a), where n = 1/2 , we show the weakly inter- 
acting, weakly perturbed case U = I^wq = 1: All MB As 
perform well. In panel b) both U and w are increased, 
but still n = 1/2: the HFA performs poorly whilst the 
other MBAs are very similar to each other and closer to 
the exact density. In panel c), n = 1/6, the HFA and 
TMA are virtually indistinguishable from the exact den- 
sity, while the BA and GWA start to deviate after some 
time. In panel d), where U = 4, = 1, none of the MBAs 
performs well. We attribute this to the poor description 
of the ground state spectral function in the band region 
by the MBA's: for w = 1, the band region provides the 
main response to the perturbation. On the other hand. 
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FIG. 4: (Color online) L = 2, n = 1, [/ = 1. a): Densities; 
exact (black) and damped GWA (thick green) densities for 
Wo = 5. b): Time reversed density, c): A(T, e) for it;o = 2 at 
T = 10 (black), T = 35 (thick green), T = 100 (thick blue), 
d): A(T,t) for = 2 at T = 10 (black), T = 35 (thick 
green), T = 100 (thick blue). 



in panel e), the TMA performs very well as the influ- 
ence of the satellite becomes essential (the other MBAs 
either lack or completely misplace the satellite, see Fig. 
|2|. As a general comment, for short ranged interactions, 
as in our system(s), the TMA is rather suitable, espe- 
cially at low densities. For the ground state, our results 
are consistent with earlier studies [25, 26 . However, the 
novel and significant aspect here is that the performance 
in equilibrium of the MBAs has considerable impact on 
the TD behavior. In particular, an external field redis- 
tributes the electrons to the unoccupied energy levels, 
and two spectral features, the energy gap and the satel- 
lite, play a key role in the TD evolution. 
MBAs and TDDFT. Our results also provide insight for 
TDDFT approaches to SCS (e.g., to the Hubbard model 
[22]). In Fig. 3, panel f), we show the KS potentials ob- 
tained via reverse engineering from the MBAs and exact 
TD densities. The performance of the different MBAs for 
Vks is consistent with the results for the TD densities. 
Time dependence U: Damping. When we prop- 
agate the KBE in our approximate, self-consistent 
schemes we obtain damped solutions and strong nu- 
merical evidence of an (artificial) steady state (Fig. 
|4^). The damped dynamics fulfills time-reversal sym- 
metry (Fig. [ZJ)). The damping increases with 
the strength of the external TD field, and is ab- 
sent in the linear response limit [27 . We investi- 
gated the instantaneous spectral function A (T, uo) = 

-Tr Im /^2T e'"^^ [G> - G<] (T + T - §) and its 
counterpart in time space, where T = (ti +t2)/2 and 
r = (ti — ^2). At the steady state, A gets broadened in 
energy space (Fig. ^) and damped in time space (Fig. 
|4]i). The time required to reach the steady state depends 
on the MBA used - in general the TMA is slowest - and 
it generally increases with system size. The damping is 
largest at the perturbed site. The steady states also de- 
pends on how the perturbation is switched on. With a 
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sudden switch-on we reach a steady state with the same 
energy as the initial (t = 0+) state. With an increasingly 
slow switch-on, the steady state gets progressively closer 
to the ground-state with a static perturbation, in line 
with the adiabatic theorem. We note that each steady 
state in a TD scheme corresponds to a solution of the 
steady KBE. Thus, in our MBAs, the steady KBE have 
multiple solutions. Some of these are artificial, they do 
not correspond to thermal equilibrium and most proba- 
bly have continuous spectral functions. 

Artificial damping in finite systems does not occur in 
wavefunct ion-based treatments. However, self consistent 
MBAs are entirely defined by a generating functional 
^[G] [1] and a corresponding set of self-consistent KBE: 
there is no direct connection to wavefunct ions. In finite 
systems, MBAs involving partial infinite-order summa- 
tions include diagrams with which annihilate more par- 
ticles than present in the system. In an "exact" theory, 
these unphysical contributions will give no net contribu- 
tion, but such perfect cancellation is in general lost in 
approximations like BA, GWA or TMA. Our numerical 
results show that these MBAs give ground state spec- 
tral functions with an incorrect pole structure, with dis- 
crete but infinitely many poles. In the time dynamics, 
these unphysical aspects of our MBAs have more drastic 
consequences. To be specific, we discuss three particle- 
conserving versions of the BA with increasing level of self- 
consistency [28]. If we evaluate the polarization P with 
ground-state propagators (Fig. [4^ top), the solution does 
not damp. If we instead evaluate P with propagators in 
the TD HFA (Fig. |4^ middle), the solution is partially 
damped. Finally, with a fully self-consistent P (Fig. |4^ 
bottom), damping occurs. Thus, if all Gs in S vary con- 
sistently with the external field, we have the remarkable 
result that a finite system reaches an (artificial) steady 
state [29,. 

In conclusion, we studied the non equilibrium dynam- 
ics of clusters with strong correlations, by comparing ex- 
act results to those from four many-body approximations 
(MBAs). Among such MBAs, the T-matrix approxima- 
tion performs very well at low densities and is in general 
superior to the GW and 2nd Born schemes. Our re- 
sults also show the principle limitations of conventional 
MBAs for finite systems, since we obtained damped solu- 
tions and artificial steady states. We attribute this to the 
inconsistency of applying infinite-order but approximate 
summations to finite systems. We acknowledge Robert 
van Leeuwen for fruitful discussions and Ulf von Barth 
for critical reading. This work was supported by the EU 
6th framework Network of Excellence NANOQUANTA 
(NMP4-CT-2004-500198) and the European Theoretical 
Spectroscopy Facility (INFRA-2007-211956). 



[1] L. P. Kadanoff and G. Baym, Quantum Statistical Me- 
chanics (Benjamin, New York, 1962) 
[2] L. V. Keldysh, JETP 20, 1018 (1965) 
[3] P. Danielewicz, Ann, Physics 152, 239 (1984) 
[4] See, for example, Progress in Non equilibrium Green's 
Functions III, J. Phys. Conf. Ser. 35, edited by M. Bonitz 
and A. Filinov (2006) 
[5] N.-H. Kwong and M. Bonitz, Phys. Rev. Lett. 84, 1768 
(2000) 

[6] A. P. Jauho, in Reference 4 , p. 313 

[7] G. Stefanucci, C.-O. Almbladh, Phys. Rev. B 69, 195318 
(2004) 

[8] N. E. Dahlen, R. van Leeuwen, and A. Stan in Reference 
4 , p. 340 

[9] M. Galperin, A. Nitzan, M.A. Ratner, Phys. Rev. B 76, 
035301 (2007) 

[10] K. S. Thygesen and A. Rubio, Phys. Rev. B 77, 115333 
(2008) 

[11] P. Darancet, A. Ferretti, D. Mayou , V. Olevano Phys. 

Rev. B 75 , 075102 (2007) 
[12] P. Myhohanen, A. Stan, G. Stefanucci and R. van 

Leeuwen, Eur. Phys. Lett. 84, 67001 (2008) 
[13] J.K. Freericks, Phys. Rev. B 77, 075109 (2008) 
[14] L. Hedin, Phys. Rev. 139, A796 (1965) 
[15] V. Galitzkii, Soviet Phys. JETP 7, 104 (1958) 
[16] D. Semkat, D. Kremp, M. Bonitz, Contributions to 

Plasma Physics, 42, 31 (2002) 
[17] K. S. Thygesen, Phys. Rev. Lett. 100, 166804 (2008) 
[18] P. Romaniello, D. SangalU, J. A. Berger, F. Sottile, L. G. 

Molinari, L. Reining and G. Onida, J. Chem. Phys. 130, 

044108, (2009) 

[19] W. Nelson, P. Pokes, Patrick Rinke and R.W. Godby, 

Phys. Rev. A 75, 032505 (2007) 
[20] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 

(1984). 

[21] R. van Leeuwen, N. E. Dahlen, G. Stefanucci, C. O. Alm- 
bladh, and U. von Barth, Time- Dependent Density Func- 
tional Theory (Springer, New York, 2006) and references 
therein 

[22] C. Verdozzi, |arXiv:0707.2317 and Phys. Rev. Lett. 101, 
166401 (2008) 

[23] H. S. Kohler, N. H. Kwong, and H. A. Yousif Comp. 

Phys. Comm. 123, 123 (1999) 
[24] C.-O Almbladh, in Reference 4 , p. 127 
[25] M. Cini and C. Verdozzi, Nuovo Cimento D 9, 1 (1987) 
[26] C. Verdozzi, R. W. Godby, S. Holloway, Phy Rev. Lett. 

74, 2327 (1995) 
[27] In this limit the response is described by the Bethe- 

Salpeter equation, with a kernel 5T^/5G. The latter would 

have a discrete spectrum in our MBAs, and so would the 

resulting density response. 
[28] U. von Barth, unpublished 

[29] This conclusion applies to all MBAs here, although par- 
tially self-consistent, but particle-conserving schemes can 
be easily devised for GWA but not for TMA. 



